This document walk you through single cell RNA-Seq analysis using various R Packages

We used sample files from this study: GSE256291. This data has not been published yet and it is used for visualization purpose only.

The data was processed using 10x Genomics protocol.

###Loading Libraries
library(rmarkdown)
library(remotes)
library(Seurat)
## Loading required package: SeuratObject
## Loading required package: sp
## 
## Attaching package: 'SeuratObject'
## The following objects are masked from 'package:base':
## 
##     intersect, t
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(presto)
## Loading required package: Rcpp
## Loading required package: data.table
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
library(devtools)
## Loading required package: usethis
## 
## Attaching package: 'usethis'
## The following object is masked from 'package:remotes':
## 
##     git_credentials
## 
## Attaching package: 'devtools'
## The following objects are masked from 'package:remotes':
## 
##     dev_package_deps, install_bioc, install_bitbucket, install_cran,
##     install_deps, install_dev, install_git, install_github,
##     install_gitlab, install_local, install_svn, install_url,
##     install_version, update_packages
library(patchwork)
library(ggplot2)

Loading datasets

First we need to load the data sets and create our Seurat object, which then let you to have all samples in one file that you can apply more downstream analysis. When we say, the data was processed using 10x technology, we have to expect 3 files for each samples, barcodes, features and matrix that has to be in one folder in order to make Read10X() function to work, otherwise they should be loaded by individualy. For more information, you could check

Now we create our Seurat objects individualy and then combine them by merge() function

seurat1 <- CreateSeuratObject(patinet1_counts, project="cd14_patient1")
seurat2 <- CreateSeuratObject(patinet2_counts, project="cd14_patient2")
seurat3 <- CreateSeuratObject(healthy1_counts, project="cd14_healthy1")
seurat4 <- CreateSeuratObject(healthy2_counts, project="cd14_healthy2")
cd14_combined <- merge(seurat1, y=c(seurat2,seurat3,seurat4), add.cell.ids = c("patient1", "patient2","healthy1","healthy2"), project = "cd14")

Once we have the objects, we can visualize the data distibution of cells in each sample. First, we calculate mitochondrian RNA to get some idea that how much of them are mitochondrian and then plot them all together. we can plot either with dots or without

cd14_combined[["percent.mt"]] <- PercentageFeatureSet(cd14_combined, pattern = "^MT[-\\.]")
VlnPlot(cd14_combined, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.

VlnPlot(cd14_combined, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size=0)
## Warning: Default search for "data" layer in "RNA" assay yielded no results;
## utilizing "counts" layer instead.

Now we want to check whether number of transcrpts and genes are correlating with each other. Alos we want to check whether mitochondrian feateares with number of RNA. This is what wer expect, positive correlation between number of features and RNAs, negative correlation between mitochondrian features and number of RNA.

plot1 <- FeatureScatter(cd14_combined, feature1 = "nCount_RNA", feature2 = "percent.mt",pt.size=1)
plot2 <- FeatureScatter(cd14_combined, feature1 = "nCount_RNA", feature2 = "nFeature_RNA",pt.size=1)
plot1 + plot2

Now we know how to filter out cells in our object, according to the QC plots that we initially generated. After that we apply normalization to our data, which is a common thing to do as we want to have all values to be comparable easier. The method of normalization is “LogNormalize” that normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor which is 10000 by default.

cd14_combined <- subset(cd14_combined, subset = nFeature_RNA > 500 & nFeature_RNA < 7000 & percent.mt < 5)
cd14_combined <- NormalizeData(cd14_combined)
## Normalizing layer: counts.cd14_patient1
## Normalizing layer: counts.cd14_patient2
## Normalizing layer: counts.cd14_healthy1
## Normalizing layer: counts.cd14_healthy2

Now we can find the most variable features, in our case genes, within our data. We can set the number of variable features as we desire, here we choose 2000. There is no a “good” number to assing here, usually between 2000 to 5000 would be informative but it depends on data and question that we have.

cd14_combined <- FindVariableFeatures(cd14_combined, nfeatures = 2000)
## Finding variable features for layer counts.cd14_patient1
## Finding variable features for layer counts.cd14_patient2
## Finding variable features for layer counts.cd14_healthy1
## Finding variable features for layer counts.cd14_healthy2

Now we can visualize the result

top_features <- head(VariableFeatures(cd14_combined), 20)
plot1 <- VariableFeaturePlot(cd14_combined)
plot2 <- LabelPoints(plot = plot1, points = top_features, repel = TRUE)
## When using repel, set xnudge and ynudge to 0 for optimal results
plot1 + plot2
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
## log-10 transformation introduced infinite values.
## Warning: ggrepel: 8 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

Since different genes have different base expression levels and distributions, the contribution of each gene to the analysis is different if no data transformation is performed. This is something we do not want as we don’t want our analysis to only depend on genes that are highly expressed. Therefore a scaling is applied to the data using the selected features, just like one usually does in any data science field.

cd14_combined <- ScaleData(cd14_combined)
## Centering and scaling data matrix
cd14_combined <- ScaleData(cd14_combined, vars.to.regress = c("nFeature_RNA", "percent.mt"))
## Regressing out nFeature_RNA, percent.mt
## Centering and scaling data matrix

Now we can apply dimensional reduction method such as PCA, UMAP and t-sne to the data so we can visualize them in 2D dimension. First we apply PCA, by default the number of top PCs is 50.

cd14_combined <- RunPCA(cd14_combined, npcs = 50)
## PC_ 1 
## Positive:  SKAP1, ETS1, RHOH, LTB, RPS27, RPS18, RPS19, RPL10, BACH2, RPS23 
##     RPS2, RPL3, RPL19, RPS12, RPS5, RPS4X, RPL10A, IKZF3, CARD11, RPL11 
##     RPL26, RPS3, SLC38A1, RPL7A, RPS6, RPL15, RPSA, RPL5, TC2N, CD247 
## Negative:  ARHGAP24, PLCB1, TRPS1, FRY, CD36, SLC2A3, TMTC2, ACSL1, PELI1, DYSF 
##     MAML3, ST6GALNAC3, SLC8A1, PADI4, S100A8, LPAR1, CYP1B1, KYNU, DIRC2, S100Z 
##     IRS2, ZNF438, GAS7, BMP2K, LIN7A, IPCEF1, JAZF1, CR1, SGMS2, SNTB1 
## PC_ 2 
## Positive:  S100A9, S100A8, AIF1, CST3, RPL26, RPL19, RPL11, RPL10, RPS23, RPL7A 
##     RPS18, S100A12, RPS12, RPS19, RPS2, HLA-DRB5, RPL15, RPS27, FCER1G, NAMPT 
##     RPS3, IER2, CFD, RPL5, RPS6, LY6E, RPL10A, S100A10, EGR1, RPL3 
## Negative:  ITGA4, HDAC9, S100Z, SKAP1, BACH2, FCHSD2, ZCCHC7, SLC38A1, DOCK10, JAZF1 
##     CARD11, ETS1, ST6GAL1, AFF1, PACS1, TMEM131L, CIITA, CBLB, IKZF3, AFF3 
##     RHOH, PCED1B, SSBP2, PDE7A, ATRN, MBD5, BCL2, RTN1, ANKUB1, IQGAP2 
## PC_ 3 
## Positive:  MTSS1, HLA-DPB1, HLA-DPA1, ITGA4, UTRN, HLA-DRA, TCF7L2, HLA-DRB1, XYLT1, FCGR3A 
##     CD74, LGALS2, SPRED1, SASH1, CIITA, CDKN1C, SLC8A1, SETBP1, KLF12, PID1 
##     HLA-DQA1, PCCA, MYOF, HLA-DMA, AIF1, AC009093.2, ST6GAL1, IFITM3, PRR5L, MGLL 
## Negative:  SLC2A3, S100A8, ACSL1, SKAP1, DYSF, PADI4, CYP1B1, S100A12, S100A9, ETS1 
##     CD247, RHOH, NAMPT, BCL11B, ITK, TC2N, IRS2, PRKCQ, INPP4B, CARD11 
##     CD96, LCK, IKZF3, CAMK4, LINC00861, ARHGAP24, PLAUR, SCML4, IL7R, BACH2 
## PC_ 4 
## Positive:  NAMPT, FCGR3A, APOBEC3A, IFI44L, CDKN1C, LY6E, EPSTI1, IFITM3, CD83, GCH1 
##     UTRN, RHOC, SETBP1, TCF7L2, OAS3, MX1, HES4, SNED1, IFI6, HERC5 
##     IFI44, MTSS1, SEMA6B, SNX9, SAT1, NFKBIZ, IFIT3, G0S2, INSIG1, PDE4B 
## Negative:  RPS4X, S100A9, RPS2, RPS6, RPL3, RPS12, S100A8, LGALS2, RPL15, RPL10 
##     RPS23, RPS18, RPL11, RPLP0, RPL19, RPS3, RPS27, RPL5, RPL26, S100A10 
##     RPL7A, RPL10A, RPSA, S100A12, EEF1B2, RPS5, CST3, RPS19, FMN1, CD36 
## PC_ 5 
## Positive:  BANK1, MS4A1, OSBPL10, EBF1, CD79A, PAX5, BLK, IGHM, SEL1L3, CD22 
##     FCRL1, TPD52, FAM129C, LINC00926, COBLL1, ANGPTL1, AFF3, PLEKHG1, BLNK, STAP1 
##     IGKC, BCL7A, FCRL2, RAB30, SPIB, TNFRSF13C, FCRL3, AC120193.1, DENND5B, ADAM28 
## Negative:  CD247, BCL11B, ITK, CAMK4, PRKCQ, LINC00861, INPP4B, IL7R, TC2N, IL32 
##     CD3D, THEMIS, CD6, CD3E, LCK, CD96, SYNE2, LEF1, NELL2, TRAT1 
##     SLFN12L, TRAC, RORA, RASGRF2, SKAP1, STAT4, TSHZ2, CD226, TXK, TRBC1
ElbowPlot(cd14_combined, ndims = ncol(Embeddings(cd14_combined, "pca")))

Here we can visualize which gene is contributing to which PC (top 4)

PCHeatmap(cd14_combined, dims = 1:5, cells = 500, balanced = TRUE, ncol = 5)

Another method of visualization of the PCA analysis.

VizDimLoadings(cd14_combined, dims = 1:2, reduction = "pca")

Another one

DimPlot(cd14_combined, reduction = "pca") + NoLegend()

cd14_combined <- RunTSNE(cd14_combined, dims = 1:20)
cd14_combined <- RunUMAP(cd14_combined, dims = 1:20)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 19:25:53 UMAP embedding parameters a = 0.9922 b = 1.112
## 19:25:53 Read 5463 rows and found 20 numeric columns
## 19:25:53 Using Annoy for neighbor search, n_neighbors = 30
## 19:25:53 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 19:25:53 Writing NN index file to temp file /tmp/Rtmp0KitSb/file170c36d0c8da5
## 19:25:53 Searching Annoy index using 1 thread, search_k = 3000
## 19:25:54 Annoy recall = 100%
## 19:25:54 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 19:25:55 Initializing from normalized Laplacian + noise (using RSpectra)
## 19:25:55 Commencing optimization for 500 epochs, with 224378 positive edges
## 19:25:58 Optimization finished
plot1 <- TSNEPlot(cd14_combined)
plot2 <- UMAPPlot(cd14_combined)
plot1 + plot2

As we see in our result, TSNE works better than the others, as we see distinct cluster among samples. Now we cherry pick 5 highly variable genes that we found in previous steps

top_features
##  [1] "FCGR3A"   "CCDC26"   "HLA-DQA1" "AFF3"     "BACH2"    "CDKN1C"  
##  [7] "CCSER1"   "TCF7L2"   "NEGR1"    "ADGRB3"   "C1QA"     "FMNL2"   
## [13] "SETBP1"   "KLF12"    "MTSS1"    "CCL5"     "NRG1"     "SNED1"   
## [19] "KCNQ5"    "HLA-DPA1"
plot1 <- FeaturePlot(cd14_combined, c("FCGR3A","CCDC26","ETS1","ITK","SEL1L3"),
                     ncol=5, reduction = "tsne")

plot1

Next step is clustering, as we need to know how many groups we can divide our samples to. A k-nearest neighbor network of cells is generated. Every cells is firstly connected to cells with the shortest distances, based on their corresponding PC values. Only cell pairs which are neighbors of each other are considered as connected. Proportion of shared neighbors between every cell pairs is then calculated and used to describe the strength of the connection between two cells. Weak connections are trimmed. This gives the resulted Shared Nearest Neighbor (SNN) network.

cd14_combined <- FindNeighbors(cd14_combined)
## Computing nearest neighbor graph
## Computing SNN

With the network constructed, the louvain community identification algorithm is applied to the netowkr to look for communities in the network, i.e. cell groups that cells in the same group tend to connect with each other, while connections between cells in different groups are sparse.

cd14_combined <- FindClusters(cd14_combined, resolution = 1)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5463
## Number of edges: 171639
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7488
## Number of communities: 14
## Elapsed time: 0 seconds
plot1 <- DimPlot(cd14_combined, reduction = "tsne", label = TRUE)
plot2 <- DimPlot(cd14_combined, reduction = "umap", label = TRUE)
plot1 + plot2

the package called “harmony”, apply various method to remove batch correction. Sometimes we need to apply it to our data when we have batch correction.

library(harmony)

cd14_combined <- RunHarmony(cd14_combined, group.by.vars = "orig.ident", dims.use = 1:20, max.iter.harmony = 50)
## Warning: Warning: The parameter max.iter.harmony is replaced with parameter max_iter. It will be ignored for this function call and please use parameter max_iter in future function calls.
## This warning is displayed once per session.
## Transposing data matrix
## Initializing state using k-means centroids initialization
## Harmony 1/10
## Harmony 2/10
## Harmony 3/10
## Harmony 4/10
## Harmony converged after 4 iterations
cd14_combined <- RunUMAP(cd14_combined, reduction = "harmony", dims = 1:20)
## 19:26:02 UMAP embedding parameters a = 0.9922 b = 1.112
## 19:26:02 Read 5463 rows and found 20 numeric columns
## 19:26:02 Using Annoy for neighbor search, n_neighbors = 30
## 19:26:02 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 19:26:02 Writing NN index file to temp file /tmp/Rtmp0KitSb/file170c364b847e5
## 19:26:02 Searching Annoy index using 1 thread, search_k = 3000
## 19:26:03 Annoy recall = 100%
## 19:26:03 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 19:26:03 Initializing from normalized Laplacian + noise (using RSpectra)
## 19:26:03 Commencing optimization for 500 epochs, with 230298 positive edges
## 19:26:07 Optimization finished
cd14_combined <- FindNeighbors(cd14_combined, reduction = "harmony", dims = 1:20) %>% FindClusters(resolution = 0.6)
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5463
## Number of edges: 181627
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7494
## Number of communities: 11
## Elapsed time: 0 seconds
DimPlot(cd14_combined, reduction = "harmony", label = F)

Next, we find top marker in each cluster and filter it the way we desire.

cd14_combined <- JoinLayers(cd14_combined)
cd14_markers <- FindAllMarkers(cd14_combined, only.pos = TRUE, min.pct = 0.25, logfc.threshold = log(1.2))
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
## Calculating cluster 6
## Calculating cluster 7
## Calculating cluster 8
## Calculating cluster 9
## Calculating cluster 10
top_marker <- cd14_markers %>% group_by(cluster) %>% top_n(n = 1, wt = avg_log2FC)
top_marker
## # A tibble: 11 × 7
## # Groups:   cluster [11]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene     
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>    
##  1 6.12e-175       1.33 0.915 0.672 2.05e-170 0       S100A12  
##  2 8.59e-123       1.60 0.618 0.277 2.88e-118 1       DYSF     
##  3 1.08e- 29       1.19 0.285 0.145 3.61e- 25 2       PALLD    
##  4 7.90e- 49       1.53 0.378 0.155 2.65e- 44 3       HLA-DQA1 
##  5 2.40e-268       4.80 0.409 0.022 8.06e-264 4       IFIT1    
##  6 3.81e- 62       3.09 0.519 0.23  1.28e- 57 5       VCAN-AS1 
##  7 0               6.22 0.561 0.017 0         6       FCER1A   
##  8 1.52e-283       6.45 0.314 0.003 5.10e-279 7       LINC02085
##  9 0              13.6  0.879 0.001 0         8       BCL11B   
## 10 0              16.5  0.333 0     0         9       IGLC3    
## 11 0              16.3  0.6   0     0         10      GP9

Now we can visualize the result

DoHeatmap(cd14_combined, features = top_marker$gene) + NoLegend()

Now we got to the most challenging part of the Single Cell analysis, cell annotation! Here, we manually annotated the cells based on marker of each cluster that we found in previous step and check the cell type based on that on panglaodb.se website, which is a great source for single cell RNA-Seq data.

new_cluster_names <- c("Macrophages",
                        "Fibroblasts",
                        "Macrophages",
                        "T Cells",
                        "Keratinocytes",
                        "Monocytes",
                        "Fibroblasts",
                        "Dendritic cells",
                        "Unknown",
                        "Epithelial Cells",
                        "T Cells",
                        "Fibroblasts",
                        "T Cells",
                        "B Cells")

names(new_cluster_names) <- levels(cd14_combined)
cd14_combined <- RenameIdents(cd14_combined, new_cluster_names)
## Warning: Cannot find identity NA
## Warning: Cannot find identity NA
## Warning: Cannot find identity NA
DimPlot(cd14_combined, reduction = "tsne", label = TRUE)

Sometimes we find to find pattern in our data, NMF clustering is a great method which allows us to find pattern from a huge data. After that we can reduce the dimension even more. We will NMF on both UMAP and TSNE

library(GeneNMF)
cd14_combined <- runNMF(cd14_combined, k = 13, assay="RNA")

cd14_combined <- RunTSNE(cd14_combined, reduction = "NMF", dims=1:13, reduction.name = "NMF_Tsne", reduction.key = "nmfTsne_")
cd14_combined <- RunUMAP(cd14_combined, reduction = "NMF", dims=1:13, reduction.name = "NMF_UMAP", reduction.key = "nmfUMAP_")
## 19:26:26 UMAP embedding parameters a = 0.9922 b = 1.112
## 19:26:26 Read 5463 rows and found 13 numeric columns
## 19:26:26 Using Annoy for neighbor search, n_neighbors = 30
## 19:26:26 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 19:26:26 Writing NN index file to temp file /tmp/Rtmp0KitSb/file170c350ff592f
## 19:26:26 Searching Annoy index using 1 thread, search_k = 3000
## 19:26:27 Annoy recall = 100%
## 19:26:28 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 19:26:28 Initializing from normalized Laplacian + noise (using RSpectra)
## 19:26:28 Commencing optimization for 500 epochs, with 211408 positive edges
## 19:26:31 Optimization finished

Now we can visualize the result from both UMAP and TSNE

DimPlot(cd14_combined, reduction = "NMF_Tsne", label=F) + theme(aspect.ratio = 1,
                                                            axis.text = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank()) + ggtitle("NMF TSNE")

DimPlot(cd14_combined, reduction = "NMF_UMAP", label=F) + theme(aspect.ratio = 1,
                                                            axis.text = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank()) + ggtitle("NMF UMAP")

Then we can see how many “real” cluster we can get from our data and then investigate the most variable gene within each cluster. According to NMF, we can have 10 cluster unline seurat clustering which was telling that we have 13 clusters.

cd14_list <- SplitObject(cd14_combined, split.by = "orig.ident")
cd14_multi <- multiNMF(cd14_list, assay="RNA", slot="data", k=4:9, nfeatures = 1000)
cd14_meta <- getMetaPrograms(cd14_multi,
                                        nMP=10,
                                        weight.explained = 0.7,
                                        max.genes=100)
NMF_Heatmap <- plotMetaPrograms(cd14_meta)

Now we zoom in each component (cluster)

lapply(cd14_meta$metaprograms.genes, head)
## $MP1
## [1] "FCGR3A" "SETBP1" "CDKN1C" "MTSS1"  "RHOC"   "TCF7L2"
## 
## $MP2
## [1] "S100A12" "RPL10A"  "RPL7A"   "EEF1B2"  "RPS4X"   "RPL15"  
## 
## $MP3
## [1] "DYSF"   "SLC2A3" "PADI4"  "MCTP2"  "VNN3"   "ACSL1" 
## 
## $MP4
## [1] "IFI44"  "MX1"    "IFI44L" "HERC5"  "OAS3"   "EPSTI1"
## 
## $MP5
## [1] "HLA-DPB1" "HLA-DQA1" "HLA-DQB1" "HLA-DPA1" "CLEC10A"  "HLA-DMA" 
## 
## $MP6
## [1] "PDE4D"  "BCL2"   "BACH2"  "ETS1"   "NEGR1"  "SEL1L3"
## 
## $MP7
## [1] "PID1"  "AFF3"  "NRG1"  "CPED1" "HDAC9" "FMN1" 
## 
## $MP8
## [1] "FOSB"       "NAMPT"      "PDE4B"      "AL390957.1" "SLC2A3"    
## [6] "CHD7"      
## 
## $MP9
## [1] "WDR49"    "RFX2"     "PDE7B"    "AGPAT4"   "FAR2"     "SERPINI2"
## 
## $MP10
## [1] "KLF12"    "HLA-DQA1" "IL1B"     "IFI44L"   "EPSTI1"   "LMNA"
library(UCell)
mp.genes <- cd14_meta$metaprograms.genes
cd14_combined <- AddModuleScore_UCell(cd14_combined, features = mp.genes, assay="RNA", ncores=4, name = "")
VlnPlot(cd14_combined, features=names(mp.genes),
        pt.size = 0, ncol=6)

sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/Toronto
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] UCell_2.8.0        GeneNMF_0.6.0      harmony_1.2.0      ggplot2_3.5.1     
##  [5] patchwork_1.2.0    devtools_2.4.5     usethis_2.2.3      presto_1.0.0      
##  [9] data.table_1.15.4  Rcpp_1.0.13        dplyr_1.1.4        Seurat_5.0.0      
## [13] SeuratObject_5.0.2 sp_2.1-4           remotes_2.5.0      rmarkdown_2.27    
## 
## loaded via a namespace (and not attached):
##   [1] RcppAnnoy_0.0.22            splines_4.4.0              
##   [3] later_1.3.2                 tibble_3.2.1               
##   [5] polyclip_1.10-6             fastDummies_1.7.3          
##   [7] lifecycle_1.0.4             globals_0.16.3             
##   [9] lattice_0.22-5              MASS_7.3-60.0.1            
##  [11] SnowballC_0.7.1             magrittr_2.0.3             
##  [13] limma_3.60.3                plotly_4.10.4              
##  [15] sass_0.4.9                  jquerylib_0.1.4            
##  [17] yaml_2.3.9                  httpuv_1.6.15              
##  [19] sctransform_0.4.1           spam_2.10-0                
##  [21] sessioninfo_1.2.2           pkgbuild_1.4.4             
##  [23] spatstat.sparse_3.1-0       reticulate_1.38.0          
##  [25] cowplot_1.1.3               pbapply_1.7-2              
##  [27] RColorBrewer_1.1-3          zlibbioc_1.50.0            
##  [29] abind_1.4-5                 pkgload_1.4.0              
##  [31] GenomicRanges_1.56.1        Rtsne_0.17                 
##  [33] purrr_1.0.2                 BiocGenerics_0.50.0        
##  [35] GenomeInfoDbData_1.2.12     IRanges_2.38.1             
##  [37] S4Vectors_0.42.1            ggrepel_0.9.5              
##  [39] irlba_2.3.5.1               listenv_0.9.1              
##  [41] spatstat.utils_3.0-5        pheatmap_1.0.12            
##  [43] goftest_1.2-3               RSpectra_0.16-2            
##  [45] spatstat.random_3.3-1       fitdistrplus_1.2-1         
##  [47] parallelly_1.37.1           DelayedArray_0.30.1        
##  [49] leiden_0.4.3.1              codetools_0.2-20           
##  [51] tidyselect_1.2.1            UCSC.utils_1.0.0           
##  [53] farver_2.1.2                viridis_0.6.5              
##  [55] matrixStats_1.3.0           stats4_4.4.0               
##  [57] spatstat.explore_3.3-1      jsonlite_1.8.8             
##  [59] BiocNeighbors_1.22.0        ellipsis_0.3.2             
##  [61] progressr_0.14.0            ggridges_0.5.6             
##  [63] survival_3.6-4              tools_4.4.0                
##  [65] ica_1.0-3                   glue_1.7.0                 
##  [67] SparseArray_1.4.8           gridExtra_2.3              
##  [69] xfun_0.46                   MatrixGenerics_1.16.0      
##  [71] GenomeInfoDb_1.40.1         withr_3.0.0                
##  [73] fastmap_1.2.0               fansi_1.0.6                
##  [75] RcppML_0.3.7                digest_0.6.36              
##  [77] R6_2.5.1                    mime_0.12                  
##  [79] colorspace_2.1-0            scattermore_1.2            
##  [81] tensor_1.5                  spatstat.data_3.1-2        
##  [83] RhpcBLASctl_0.23-42         utf8_1.2.4                 
##  [85] tidyr_1.3.1                 generics_0.1.3             
##  [87] S4Arrays_1.4.1              httr_1.4.7                 
##  [89] htmlwidgets_1.6.4           uwot_0.2.2                 
##  [91] pkgconfig_2.0.3             gtable_0.3.5               
##  [93] lmtest_0.9-40               XVector_0.44.0             
##  [95] SingleCellExperiment_1.26.0 htmltools_0.5.8.1          
##  [97] profvis_0.3.8               dotCall64_1.1-1            
##  [99] Biobase_2.64.0              scales_1.3.0               
## [101] png_0.1-8                   spatstat.univar_3.0-0      
## [103] knitr_1.48                  rstudioapi_0.16.0          
## [105] reshape2_1.4.4              nlme_3.1-164               
## [107] cachem_1.1.0                zoo_1.8-12                 
## [109] stringr_1.5.1               KernSmooth_2.23-22         
## [111] parallel_4.4.0              miniUI_0.1.1.1             
## [113] pillar_1.9.0                grid_4.4.0                 
## [115] vctrs_0.6.5                 RANN_2.6.1                 
## [117] urlchecker_1.0.1            lsa_0.73.3                 
## [119] promises_1.3.0              xtable_1.8-4               
## [121] cluster_2.1.6               evaluate_0.24.0            
## [123] cli_3.6.3                   compiler_4.4.0             
## [125] crayon_1.5.3                rlang_1.1.4                
## [127] future.apply_1.11.2         labeling_0.4.3             
## [129] plyr_1.8.9                  fs_1.6.4                   
## [131] stringi_1.8.4               viridisLite_0.4.2          
## [133] deldir_2.0-4                BiocParallel_1.38.0        
## [135] munsell_0.5.1               lazyeval_0.2.2             
## [137] spatstat.geom_3.3-2         Matrix_1.7-0               
## [139] RcppHNSW_0.6.0              future_1.33.2              
## [141] statmod_1.5.0               shiny_1.8.0                
## [143] SummarizedExperiment_1.34.0 highr_0.11                 
## [145] ROCR_1.0-11                 igraph_2.0.3               
## [147] memoise_2.0.1               bslib_0.7.0